Numerical construction of a low-energy effective Hamiltonian in a self-consistent 
Bogoliubov-de Gennes approach of superconductivity 
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We propose a fast and efficient approach for solving the Bogoliubov-de Gennes (BdG) equations in 
superconductivity, with a numerical matrix-size reduction procedure proposed by Sakurai and Sug- 
iura [J. Comput. Appl. Math. 159, 119 (2003)]. The resultant small-size Hamiltonian contains the 
information of the original BdG Hamiltonian in a given energy domain. In other words, the present 
approach leads to a numerical construction of a low-energy effective theory in superconductivity. 
The combination with the polynomial expansion method allows a self-consistent calculation of the 
BdG equations. Through numerical calculations of quasi-particle excitations in a vortex lattice, 
thermal conductivity, and nuclear magnetic relaxation rate, we show that our approach is suitable 
for evaluating physical quantities in a large-size superconductor and a nano-scale superconducting 
device, with the mean-field superconducting theory. 

PACS numbers: 74.20.Rp, 74.25.0p, 74.81.-g 



I. INTRODUCTION 

To solve an eigenvalue equation is one of the central 
issues in condensed matter physics. The ground state in 
a many-body system is nothing but the eigenvector asso- 
ciated with the lowest eigenvalue of a many-body Hamil- 
tonian. The Lanczos algorithm in the exact diagonaliza- 
tion ^ is suitable for this issue. The critical temperature 
in superconductivity is evaluated by the greatest eigen- 
value of the linearized Eliashberg equations. The power 
iteration algorithm is useful for solving these equations. 
Thus, a lot of efhcient methods for either minimum or 
maximum eigenvalues have been developed. 

In superconductors, low-energy quasiparticle excita- 
tions are quite important for examining thermodynamic 
quantities, transport properties, and so on. Their en- 
ergy scale is characterized by a superconducting gap en- 
ergy (~ meV), much smaller than a band width (~ eV). 
In the mean-field Bardeen-Cooper-SchriefFer (BCS) the- 
ory, these excitations correspond to the eigenvalues at 
the center of an energy distribution of the Bogoliubov- 
de Gennes (BdG) Hamiltonian^.. This is a direct conse- 
quence of the particle-hole symmetry of the BdG Hamil- 
tonian. Furthermore, such an intermediate region can 
include zero eigenvalues (i.e., zero modes). The zero 
modes are related to fundamental properties of topolog- 
ical insulators and superconductors^. Therefore, in or- 
der to study bulk properties of various superconductors 
and nano-scale devices with topological materials from 
atomic-scale physics, an efficient method to obtain an 
intermediate spectral region of the BdG Hamiltonian is 
highly desirable. 

Typically, the full diagonalization method is used for 
solving the BdG equations—"—. However, this approach 
requires a lot of computational memories and a long com- 



putational time. In contrast, the polynomial expansion 
method—"— allows efficient self-consistent calculations in 
superconductivity, without any diagonalization. This ap- 
proach drastically reduces a computational cost and has 
an excellent parallel efficiency. However, it is not easy to 
directly obtain the eigenvalues and the eigenvectors of the 
BdG Hamiltonian. Thus, this method is not suitable for 
calculating dynamical correlation functions (two-particle 
Green's functions), which lead to important quantities 
such as spin/charge susceptibilities, nuclear magnetic re- 
laxation rate, optical/thermal conductivities. Hence, the 
algorithms for treating an intermediate energy region of 
the BdG Hamiltonian have not been adequately studied. 

In this paper, we propose a fast and efficient method 
for numerically calculating the eigenvalues and the eigen- 
vectors of the BdG equations. Our approach is the 
combination of the polynomial expansion method with 
a contour-integral-based method developed by one of 
the present authors (TS) and Sugiura (Sakurai-Sugiura 
method)^^"— . The Sakurai-Sugiura (SS) method allows a 
systematic reduction of a matrix size of a target Hamilto- 
nian, with keeping the information relevant to low-energy 
excitations in a superconductor. A contour-integral rep- 
resentation of the projection operator onto an energy do- 
main plays a crucial role. We obtain an effective Hamil- 
tonian only with a gapless surface state in topological in- 
sulators and superconductors in large-scale systems, for 
example. 

Let us summarize our approach for calculating phys- 
ical quantities in the mean-field superconducting the- 
ory. First, we perform a self-consistent calculation of the 
BdG equations to obtain a superconducting gap func- 
tion. Next, we numerically derive a low-energy (small- 
size) effective Hamiltonian from the BdG Hamiltonian 
with the resultant superconducting gap. Finally, we cal- 
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culate physical quantities using the eigenvalues and the 
eigenvectors of this effective Hamiltonian. 

This paper is organized as follows. In Sec. HH we show 
a general formulation of a mean-field fermionic theory. 
The Green's functions are expressed by the eigenvalues 
and the eigenvectors of the BdG equations. In Sec. IIIIl 
we explain the polynomial expansion scheme. The appli- 
cation of the SS method to superconductivity is proposed 
in Sec. IIVI We show the theoretical background of this 
approach and the algorithm. We stress that the present 
algorithm is suitable for parallel computation since the 
procedure is composed of solving a set of the linear equa- 
tions which are independent of each other. In Sec. |Vl 
we show the results for typical examples, as well as the 
computational costs of the present proposal. We perform 
large-scale calculations in various physical situations. As 
for inhomogeneous systems, we consider a vortex lat- 
tice on a two-dimensional square lattice. The quasiparti- 
cle excitation spectrum is obtained, varying a magnetic 
field and a coherence length. The thermal conductivity 
is also evaluated. Moreover, we examine temperature- 
dependence of the nuclear magnetic relaxation rate in an 
s-wave superconductor, as an example of a uniform su- 
perconductor. These demonstrations indicate that the 
present approach is a fast and accurate method for nu- 
merically constructing a low-energy effective Hamiltonian 
in the mean-field superconducting theory. Section IVII is 
devoted to the summary. 



II. FORMULATION 
A. Hamiltonian 

Throughout this paper, we set ft = /cb = 1. Let us con- 
sider a Hamiltonian for a fermionic many-body system, 

with c = (ci,C2, . . . ,CAr)'^ and c = (c| , cj, . . . , cjy)'^. 
Here, the symbol T represents transposition. The 
fermionic annihilation and creation operators are denoted 
as, respectively, Ci and c| {i = 1,...,7V). The index i 
includes all the relevant degrees of freedom such as spa- 
tial sites, spins, orbitals, and so on. The canonical anti- 
commutation relation is [ci,cj]+ — 6ij. The Hamiltonian 

matrix _ff is a 2N x 2N Hermite matrix. The hermitian 
property of H and the canonical anti-commutation re- 
lation imply that the N x N complex matrices A and 
B in H satisfy A' = A and ^ —B. In the case of 
superconductivity, H corresponds to the mean-field BCS 
Hamiltonian and B contains superconducting gaps. 



B. Bogoliubov-de Gennes equations 

The BdG equations are regarded as the eigenvalue 
equations with respect to H 

Hx^ = e^x^, (7 = 1,2,--- ,2Af), (2) 

with — (it-y, D-y)"'". The column vectors and 
are iV-component complex vectors. To solve the BdG 
equations is equivalent to the diagonalization of H with 
a unitary matrix U . The matrix elements of U are 

C^i,7 = 1*7, ij Ui-\-N^^ = Vj,i. (3) 

The eigenvalues are not independent of each 
other. In fact, using the particle-hole transforma- 
tion— such that Jx^ — {v*, u*)"^ , one can show that 

H{JXj) — —ej{Jxj). 

C. Two- particle Green's functions 

Various important physical quantities in condensed 
matter physics are related to two-particle Green's func- 
tions. Here, we write them in terms of the solutions of the 
BdG equations. Let us consider a two-particle Green's 
function in the imaginary time r, 

Qi234(t) - (T,[c|^(t)c,,(t)c1^(0)c.,(0)]), (4a) 
= G,,,3(t)G,,,,(t) - F,,,Jr)i^,,,3(T), (4b) 

with the one-particle Green's functions 

G,,(T) = -(T,[c,(r)c](0)]), (5a) 
F,,{T) = -{Tr[c,iT)c,m), (5b) 

F. ,{T) = -{TA4{r)c]m), (5c) 

G, j(t) = -(T,[4(t)c,(0)]). (5d) 

With the use of the relation 

/ dre^"-^A(r)S(r) = ^ V A(^c^„)S(^^!„, - zc^„), 

(6) 

the Fourier transformed function Q(if2m) is 

-Fi^i^{iuJn)Fi^iS^m - il^n)] • (7) 

Here, /3 is the inverse temperature and a;„ = 7r(2n-|- 1)//3 
and = Tr(2m)/f3 are the fermionic and bosonic 

Matsubara frequencies, respectively. The one-particle 
Green's functions are written as a 2N x 2A^ matrix, 

(8) 
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We find that = , Fij = G^.j+n, Fij = Gi+Nj, and 
Gij = Gi+N,j+N ■ We can phenomenologically describe a 
dissipation efi^ect, replacing the delta function in A with 
an approximate J-function. The dynamical correlation 
function with the real energy Q. is 



2^ 



Ql234(f^) — 2J ^»2,7^«i+JV,7' \Ui3,-tUii+N,^' 
1,1' 



with setting i^m ^ + if] {rj ^ 0+). Here, f{x) 
l/{e^^ + 1) denotes the fermion distribution function. 



III. POLYNOMIAL EXPANSION METHOD 

We briefly summarize the polynomial expansion 
method for a self-consistent calculation of the BdG equa- 
tions, according to our previous paper^^. The essence is 
the expansion of the spectral density of the Green's func- 
tions, with orthonormal polynomials in [—1, 1] satisfying 



5[X - x') = ^^^^ (j)n{x)(l)n {x'), 

^ — ^ ?/i _ 



n=0 
1 



W„(S„,m = / (l)n{x)(l)m(x)W(x)dx, 



(10a) 



(10b) 



4>n+l{x) = (a„ + hnX)(t>n{x) - C„(/)„_i (x). (lOc) 

The Chebyshev polynomial is often used, because the 
resultant formulae become the simplest forms. The ap- 
plication of the other polynomials is discussed in, e.g., 
our previous paper . 

The spectral density (matrix) is given as a difference 
between the retarded and the advanced Green's func- 
tions, d{uj) = G{uj + iQ) — G{(jj — iG). Let us expand 
d{iS) by (t>n{x), rescaling H and a; so that K. = {H — b)/a 
and X = {uj — b)/a, with a = (E'max — £'min)/2 and 
b = (-Emax + £^min)/2- Hcrc, E'niax and £;,nin are energy 
scales satisfying i?niin < £7 < -Emax- The elements of 
d{uj) are related to various correlation functions. Using 
the constant vectors e{i) and h{i) such that [e(i)]^ = Si,^ 



and [h{i)]j — Si+N,-y, we obtain 



where 



00 ^ 



r„=/ dxf{ax + b)W{x)(j3n{x). 



(11a) 
(lib) 

(12) 



A sequence of a vector q„ = (j)n{IC)q [q = e{i), h{i)] is 
recursively obtained by 

Qn+l = {an + bnVjQn - c„(7„_i {n > 2), (13a) 
qi=0i(/C)q, qo^(j}Q{K)q. (13b) 

The use of the recurrence formula leads to a self- 
consistent calculation of the BdG equations, without any 
diagonalization of H . 



IV. CONTOR-INTEGRAL-BASED METHOD 
(SAKURAI-SUGIURA METHOD) 

In this paper, we use the SS method for finding eigen- 
values in a given energy domain and their associated 
eigenvectors. This approach is a numerical solver for a 
generalized eigenvalue problem so that Ax = XBx, with 
A, B E C"=^"% and has been applied to various physical 
issues such as the real-space density functional theory— 
and the lattice quantum chromodynamics^*^. In this pa- 
per, B is the identity matrix, and A is an Hermite matrix. 

Our aim is to reduce the size of A, keeping as much 
information of the eigenvalues and the eigenvectors as 
possible. Let us consider the use of an x mg {ug > rris) 
matrix Q, whose ith. column is qi G C"" (i.e., Q = 
{qi, . . . , Qto^}). Here {<7i}™\ is a set of hnearly inde- 
pendent vectors in C"' . We obtain an rus x TOs matrix 



A = Q^AQ. 



(14) 



We denote a given energy domain of A as ^r(c K). Let us 
suppose that qi is represented by a linear combination of 
{xjYJl^^, where Axj = XjXj {Xj G £). Thus, A contains 
nis eigenvalues of A ({Aj}™!^). It is necessary for imple- 
menting this procedure to know parts of the eigenvectors 
of A. Furthermore, one has to carefully choose mg to 
avoid losing the relevant information of A. Remarkably, 
this issue will be solved by an approximate evaluation of 
contour integrals associated with a projection operator 
onto eigenspaces of A. All the steps of the algorithm is 
summarized in Sec. lIVD] 



A. Projection and moment vectors 

We start with a way to make a projection onto a target 
subspace spanned by {xj}^^ iW^jW = !)■ An arbitrary 
ng-dimensional vector v is expanded by {xi}, with rtg 
complex coefficients, v = J^'^i^i- We define the projec- 
tion Pr{A) as 



PriA)v ^"^ajXj. 



(15) 



Using Pi = Xix\, we find that the resolvent-- of A is 



p,. 



zl - A z-u 
1=1 



(z e CV(A)), (16) 
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with the Us x Us identity matrix, / and a set of ah the 
eigenvalues of A, a (A), since A = J2^iPi- Let us suppose 
that the rUs distinct eigenvalues (i.e., simple poles on C) 
are located inside a closed loop F on C, and the others 
are outside F, as shown in Fig. [T] Thus, we obtain a 
contour-integral representation of Pr (A) 



PriA) 



dz 



1 



2m zl - A 



(17) 



Now, let us write essential quantities for determining 
the reduction matrix Q. The moment vector (k = 
0, 1, . . . , M - 1) is defined as Sk = A''Pr{A)v, with a 
vector V G C"" . From its contour-integral representation, 
we find that Sk is related to the fcth moment. 



dz 



r 2TTi zl - A 



(18) 



An important property of Sk is that this is a vec- 
tor in the eigen-space associated with Py{A). This 
fact is checked by applying A'' to Eq. (fTSt . We 



remark that all the moment vectors are linearly in- 
dependent of each other for an arbitrary v, since 
the subspace spanned by {sfcl^fjo^ is the order-M 
Krylov subspace^^ generated by A, K,m{AtPt{A)v) = 
span{Pr(A)i;, APt{A)v, A^Pt{A)v, • • • , A^-^Pt{A)v}. 
In our algorithm to determine Q, M is an input param- 
eter. Then, one has to construct linearly independent 
vectors from {s/cj^lo^' varying v, and evaluate mg with 
a proper manner. Another important issue is to numer- 
ically calculate the counter integrals. These points will 
be explained in the following. 



Im z 



Re z 



the moment vector is approximately written by 



1 



(19) 



with Wj ~ w{Oj), w{9) = -~iC,'{9), Zj — j + p(j, and 
Cj ~ The vector yj is the solution of the linear 

equation [zjl — A)yj = v. When all the eigenvalues are 
located on the real axis, it might be better to put the 
quadrature points closer to the real axis, 



2n 



1 



+ p(cos Oj + ia sin 9^ ) , ^ [j - 2 J ' ^^^'^ 

with vertical scaling factor a (0 < a < 1). The quadra- 
ture weight is 



Wi 



a cos 9j + I sm ( 



(21) 



When a = 1, Fq is a unit circle. Other formulae includ- 
ing contour integrals are also calculated by this A'q-point 
numerical quadrature. 



C. Construction of subspace 

Now, we show a way to determine the subspace size 
and the corresponding reduction matrix. First, we make 
a sequence of the moment vectors, varying v. The use 
of this sequence is essential for constructing the linearly 
independent vectors and the subspace. We set L complex 
vectors, z;*(g C"=) (i = 1, 2, . . . , L), and make an rig x 
L real matrix V = {v^ . ^ . . . , v^}. We call V source 
matrix. Then, we obtain an ng x L matrix 5^-, 



1 



with 



{z,I - A)Y, = V. 



(22) 



(23) 



FIG. 1. Schematic diagram of a contour on 



B. Approximation of contour integrals with 
numerical quadrature 

We show a method to approximate a contour integral 
with numerical quadrature. Let us suppose that a Jor- 
dan curve F on C is represented by scaling and shifting 
another Jordan curve Fo, with a scaling factor p and a 
shift 7. Without loss of generality, we assume that Fo 
encloses the origin on C. Let (^{9) be a point on Fq, with 
a parameter 9 {0 < 9 < 27r), and let z on F be given by 
z{9) = "f + p(^(9). Then, using A^^-point quadrature rule, 



The ith column of Sk is related to v'-, s^. = 
i^/Nci)J2jPWjz^y], with {zjl - A)y] = v\ Each el- 
ement of is a uniform random variable in (—1, 1). It 
means that we make A*'Py{A) via random sampling, with 
the source size L. The integer parameter L is determined 
by the calculation of Tr Pr {A) , as seen below. 

Next, we determine the subspace size. A bundle of Sk 
leads to an Ug x LM matrix S — {S'o, S*!, . . . , S'a/-i}. 
Now, we perform the singular-value decomposition of S, 
and obtain the singular values {cJi}, with cti > (72 > 
. . . > 0. Then, we find the number of the singular val- 
ues satisfying ai/cri > S, with a small positive constant 
10"^''). Thus, we have an effective rank (i.e., the 
number of the predominant linearly independent vectors) 
of S. We stress that this rank should be greater than or 
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equal to a prior value rhg, which is determined by cal- 
culating the trace of the resolvent matrix (see below). 
We use the resultant effective rank as the dimension of 
the subspace. We remark a large source size leads to a 
highly accurate calculation, but causes an extreme in- 
crease in the subspace size. To avoid this increase for the 
high accuracy, one can use an iterative refinementlA of a 
subspace, as seen in Appendix [Xj 

Now, the construction of the reduction matrix to the 
subspace is straightforward. Using a submatrix com- 
posed of the first mg columns of (e.g., Sk{-, 1 : mg) 
in terms of Fortran) and the Gram-Schmidt orthonor- 
malization, we obtain an Us x mg matrix Q whose mg 
column vectors are orthonormal to each other. This is 
our reduction matrix. From the construction manner, 
one finds that Q contains predominant rus eigenvectors 
of A. Alternatively, one may use a matrix U such that 
S = Utw'f and t = diag(cri,cr2, ■ • ■)• This matrix IS au- 
tomatically obtained when performing the singular-value 
decomposition of S by using the ZGESVD routine of LA- 
PACK, and the iris column- vectors are orthogonal to each 
other. 

The source size L and the moment size M have to be 
carefully chosen. In particular, the integer LM, which 
is the total number of the moment vectors to take in a 
simulation, should be as small as possible for a few com- 
putational costs. First, we predict the prior rank rhg, 
with the stochastic estimation methodi^ii^. We prepare 
an Us X Lq real matrix V whose elements are either — 1 
or 1 with equal probability. Here, Lq is an input parame- 
ter. The number of the eigenvalues inside F is Tr Pr{A). 
The stochastic estimation of TrM for an x ma- 
trix is TtM - {l/Lo)Y.f=i{v'VMv'{See, Appendix|B|. 
Thus, using the A'q-point numerical quadrature, rhs is 
estimated by 

rfig^^Y^iv-fs^o- (24) 

^0 z=l 

Then, the source size L is 



with K > 1. The symbol [x] means the smallest integer 
greater than a; . The value of LM is larger than mg. Thus, 
the requirement that the rank of S is greater than and 
equal to rhs is automatically satisfied. 

D. Algorithm of the SS method 

Now, we show all the steps of the SS method for cal- 
culating the eigenvalues and the eigenvectors of the BdG 
Hamiltonian H. 

(i) Set H e C^^^ {us = N), Lq, M, and 
V — {v^, ■ • • , v^°}. The elements of the sampling vec- 
tor Vi take either -1 or 1, with equal probability. 



(ii) Solve Eq. ^ for Y,-, j = 1, • • • , A^q. One can solve 
these equations separately so that parallel computations 
can be easily implemented. 

(iii) Compute Eq. ((22)) . 

(iv) Compute Eq. ([M)) . and estimate L via Eq. (1^ . 

(v) Give the elements oiV — {v^, • • • , v^} by random 
numbers and solve Eq. ([25)1 . 

(vi) Compute Eq. ([^^ using the results in (v). 

(vii) Perform the singular-value decomposition 
U'EW'^ = {So, ■ ■ ■ , Sm-i} and find nis such that 
Cj I / 1 ci I < 6 for 1 < j < rUs. 

(viii) Obtain a matrix Q from Q — U{:,1 : rris). 

(ix) Form H ^ Q'^HQ. 

(x) Compute the eigenvalues ti and eigenvectors Wi of 
the matrix H . 

(xi) Set Xi — Qwi. 

If one uses the iterative refinement of a subspace, one 
adopts either Eqs. (|A1I) or (jA3l) . and goes to (vi). The 
one- or two- particle Green's function are calculated by 
the eigen-pair {ei,Xi). 



V. NUMERICAL DEMONSTRATIONS 

We show the effectiveness and the validity of the 
present approach, focusing on a single-band superconduc- 
tor. Hereafter, the index i in Sec. Ill Al indicates a spatial 
site on a two-dimensional square lattice. The spin indices 
('|~ and D are explicitly written in the creation and annihi- 
lation operators. We consider a two-dimensional Lx x Ly 
lattice system, with the nearest-neighbor hopping t. The 
spatial site index runs from 1 to L^; x Lj,. We impose the 
periodic boundary condition. The superconducting gap 
equations is given as A.^ — Vij{ci i^Cj^^), with pairing in- 
teraction Vij . This equation is self-consistently solved by 
the polynomial expansion method, as shown in Sec. IIIII 
The parameters in the SS method are set as Lq ~ 10, 
M = 16 and k = 1.5. We remove the eigen-pair {ei,Xi) 
whose relative residual is greater than 10^^, as spurious 
eigenpairs. See Eq. (^5)) . We do not use the iterative 
refinement of a subspace in the following examples. 



A. Computational costs in self-consistent 
calculations 

We use the polynomial expansion scheme to obtain 
self-consistent superconducting gaps. Let us evaluate 
computational costs for a d-wave superconductor at zero 
temperature, on a square lattice {L^ = Ly). An eval- 
uation for an s-wave superconductor was shown in our 
previous contribution—. We measure the elapsed time 
for calculating and updating the order parameters at one 
iteration. 

We used a superconiputing system PRIMERGY 
BX900 in Japan Atomic Energy Agency. As shown 



6 




System size N 



FIG. 4. (Color online) Order-parameter distribution obtained 
FIG. 2. (Color online) System-size dependence of elapsed by a self-consistent calculation with the polynomial expansion 
time for calculating and updating order-parameters at one scheme for a vortex lattice in an s-wave superconductor, with 
iteration with the polynomial expansion scheme for a d-wave on-site hopping Vu — — 1.5t and chemical potential fj, = 0. 
superconductor at zero temperature, on a L^, x Ly square 
lattice {Lx = Ly). The system size is N = Lx x Ly. 

B. Eigenvalues in a vortex lattice system in an 
s-wave superconductor 
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FIG. 3. (Color online) Strong-scaling plot in self-consistent 
calculations with the polynomial expansion scheme. 



in Fig. [2J the elapsed time for one iteration is sub- 
jected to an 0{N'^) rule, with increasing the system size 
N = 2{L^ X Ly). This tendency is kept from 32 to 4096 
CPU cores. In contrast, the full diagonalization scheme 
inevitably demands 0{N^) costs in the core part of a 
calculation. This is a big advantage of the polynomial 
expansion scheme. Furthermore, we focus on the strong 
scaling, as shown in Fig. [3] One can see an excellent 
strong scaling up to 4096 CPU cores. 



We show the eigenvalues obtained by the SS method. 
The system in this section has a vortex square lattice in 
an s-wave superconductor. The parameters are set as fol- 
lows: on-site interaction Vu — —1.5t, chemical potential 
/i = 0, and spatial size Lj^ x Lj, = 64 x 64. The ma- 
trix dimension is TV = 8192 with 49152 nonzero entries. 
The rescaling parameters in the polynomial expansion 
method are a = 8t and 6 = 0. Also, a cut-off parame- 
ter in the polynomial expansion schemeil is 2000. The 
resultant order parameter is shown in Fig. |4l 

The relative residual for the eigen-pair {ei,Xi) is cal- 
culated by 



resi 



\Hx, 



Hxi 



(26) 



After the self-consistent calculations with 100 iterations, 
we obtain the eigen-pairs with the use of the SS method. 

The quadrature points are set by Eq. (l20t and the cor- 
responding weights are set by Eq. (PT|) . with a — 0.5. 
The contour F is set as 7 = and p = O.lbt. We 
adopt a sparse solver PARDISO^"^ to compute Eq. (^5]) . 
This solver uses the nested dissection algorithm from the 
METIS package^i. We also confirm that the calculations 
with the shifted BiCG method, which is suitable for large 
sparse matrices, have a similar result. 

Let us investigate iVq-dependence of the eigenvalues 
and the relative residual. Figure [S] shows that a calcula- 
tion with A^q — 64 has good precision about the eigenval- 
ues inside F (—0.15 < < 0.15). We note that the cal- 
culation with A'q = 64 takes about 30 .seconds with only 
one CPU core (Intel Xeon X5550 2.66GHz) by a desktop 
computer. The conventional full diagonalization takes 
about 40 minutes with the same machine. When the 
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(a) eigenvalues 



0.15 




5 10 15 20 25 30 35 40 45 50 
index 



FIG. 5. (Color online) (a) Eigenvalues of the Bogoliubov- 
de Gennes Hamiltonian for the order parameter shown in 
Fig. |4l changing the total number of the quadrature points 
Nq. The horizontal axis represents the index of the eigenvalue, 
in ascending order, (b) Relative residual for each eigen-pair 
(ei, Xi), varying A^'q. The index i in the horizontal axis is the 
same as in (a). 



system size becomes quadruple {Lx x Ly = 128 x 128), 
it takes about 5 minutes to obtain the same eigenvalue 
distribution, with the same one CPU core. 

We discuss the accuracy of an eigenvalue calculation 
in the SS method. This issue depends on the param- 
eters related to the contour integral representation, as 
well as the source size L. The simplest improvement can 
be achieved by increasing the total number of the quadra- 
ture points, A^q. Alternatively, we obtain better accuracy 
with smaller iVq, continuously deforming the contour T. 
Figure [6] shows the relative residual, varying TVq and the 
vertical scaling factor a. We find that the accuracy be- 
comes higher than Fig. [5jb), even though TVq is small. 
Here, we take a relatively larger source size (k = 2 and 
M = 10 in Eq. (|25|)). compared with the previous calcu- 
lations. 




index 



FIG. 6. (Color online) Accuracy of the Sakurai-Sugiura 
method in terms of relative residual H26p . The horizontal axis 
corresponds to the index of the eigen-pair, as seen in Fig. [5] 
The physical parameters are the same as in Fig. U but the 
parameters related to the numerical calculations are changed. 
The total number of the quadrature points for contour inte- 
grals is set as either Nq — 32 or 64. The source matrix size L 
in Eq. (|25|) is changed by varying a numeric constant k and 
the moment size M. The contour in the contour integrals is 
continuously deformed by the vertical scaling factor a, as seen 
in Eq. ((20)) . 



C. Magnetic-field dependence of tlie eigenvalue in 
long-coherence-length superconductors 

We show the magnetic-field dependence of the eigen- 
values in a long-coherence-length superconductor. A co- 
herence length of a superconductor f is roughly estimated 
by the ratio of the Fermi velocity to the amplitude of the 
order-parameter ^ vp/A). In many materials expect 
for high-Tc cuprates, the coherence length is much larger 
than an atomic length. Typically, the electric states in 
such superconducting systems are described by the qua- 
siclassical Eilenberger theory^^, with neglecting atomic- 
scale physics. However, interesting microscopic phenom- 
ena such as an interference effect and discretized quan- 
tum bound states in a vortex core are never treated in 
this approach. The mean-field BdG approach can treat 
these atomic-scale phenomena, but requires extremely 
large computational costs when the coherence length is 
much larger than an atomic length. Thus, to solve the 
BdG equations in a long-coherence-length superconduc- 
tor is a challenging issue. 

We consider a vortex lattice in an s-wave superconduc- 
tor, with on-site interaction Vu = —2t and chemical po- 
tential n = —t. These parameters correspond to a model 
with a relatively smaller Fermi surface, compared with 
Sec. lVBl The temperature is set as T = 0.04i. We use the 
domain F with 7 = and p = O.lt (-0.lt < < O.K), 
and iVq = 64. The magnetic field becomes small with 
decreasing the system size, since the total magnetic flux 
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FIG. 7. (Color online) Order-parameter distribution ob- 
tained by a self-consistent calculation with the polynomial 
expansion scheme for a vortex lattice in an s-wave supercon- 
ductor, with on-site hopping Vu = ~2t, chemical potential 
fi — —t, temperature T = 0.04t. The spatial size is given by 
— Ly — 64. 



is fixed. As shown in Fig. [71 the amplitude of the order- 
parameter is similar to that in the previous section. Com- 
bined with the small Fermi surface, the coherence length 
is longer than in Sec. VB. Under the periodic boundary 
condition, two bound states appear at each vortex core. 
Their energy eigenvalues are degenerate when the system 
size is large (low magnetic field). With increasing the 
magnetic field, a splitting in the degenerate eigenvalues 
occurs, as shown in Fig, [5l^a). The splitting comes from 
the occurrence of an overlap between the bound states in 
a vortex core. Figure[5]Jb) shows quantum oscillation as a 
function of the inter-vortex distance originating from an 
interference effect between two vortex bound states. We 
note that it is hard to discuss the degeneracy splitting 
with the only use of the polynomial expansion scheme, 
since the polynomial expansion calculates the local den- 
sity of the states not the eigenvalues. 



D. Magnetic field dependence of tlie thermal 
conductivity 



We show the magnetic-field-dependence of the thermal 
conductivity in an s-wave superconductor, with a vortex 
lattice. All the physical parameters are the same as the 
ones in Sec. IV CI We use the domain T with 7 = and 
p = 0.45i (-0.45i < < 0.45t), and = 64. It takes 
about one hour to obtain 1857 eigenvalues located in the 
domain F with the same one CPU core when the system 
size is Lx Ly — 90 X 90 as shown in Fig. O One can 
clearly find that the gap amplitude in the bulk states is 
Aq ~ 0.2t, which is consistent with the estimation with 
the use of Fig. H 

Let us consider the case when a temperature gradient 
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FIG. 8. (Color online) Magnetic-field dependence of the dis- 
tribution of the eigenvalues as functions of (a) the eigenvalue 
index and (b) the inter-vortex distance. 



exists along x-axis. Using the linear response theory—, 
the electronic thermal conductivity per volume of a su- 
perconductor is 



= — lim — Im \Rr 



n + io)] 



l,J2F,y\[U^Vm]y,\\ 



(27) 
(28) 



7,7' 



with 



F , — 
^77 ~ 



xn5{uj - uj')uj'^f'{uj). 



(29) 



Here, P^xii^m) is the Fourier transformation of a 
current-current correlation function— '^^1^^ in the imagi- 
nary time, Pxx{t,0) = (T^[Ju,£c(t) Ju,£c(0)]), where Ju,x is 
the Heisenberg operator of energy flux^^ along a;-axis. 
The Lorentian kernel S,f{uj) — ri/[Tr{uj'^ + rj'^)] in F-^^i 
represents a dissipation effect^. The matrix U contains 
the eigenvectors of the BdG Hamiltonian, while the ma- 
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FIG. 9. (Color online) Eigenvalue distribution in an s-wave 
superconductor with a vortex lattice. The spatial size is given 
by L^ = Ly = 90. 



FIG. 11. (Color online) Temperature dependence of the nu- 
clear magnetic relaxation rate in a uniform s-wave supercon- 
ductor. Lx = Ly = 64. 
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FIG. 10. (Color online) Magnetic-field dependence of the 
thermal conductivity as a function of an inter- vortex distance. 
The inter- vortex distance become long, when decreasing mag- 
netic field. 



trix V'^^'> includes contributions from the energy flux. We 
show their explicit formulae in Appendix [Cl 

We adopt the damping factor rj — O.OOSt. We as- 
sume that 9ij = in Eq. (jC2[) because the vector po- 
tential around a vortex is small. Figure [TU] shows that 
Kxx drastically increases when the inter-vortex distance 
is shorter than around 60 (i.e., a high magnetic- field do- 
main). This behavior could be related to an interference 
effect between the two bound states in vortex cores in 
high magnetic field, as seen Fig. [T] A quantum oscilla- 
tion in the eigenvalue distribution becomes remarkable 
in such a high magnetic-field domain. 



E. Temperature dependence of nuclear magnetic 
relaxation rates 



We show the temperature dependence of the nuclear 
magnetic relaxation rate in a uniform s-wave supercon- 
ductor. It is well known that the nuclear magnetic relax- 



Ti(n,r) = 



i?(r.,r)' 



i?(n,r)= linr Im^^^il^, 



(30) 
(31) 



Q0,eQ>O,e3>O 

X TTTf'{e^)5{e^ - ep). (32) 

We use the eigenvalues in the domain F with 7 = 0.25i, 
p = 0.25t (0 < Cq < 0.5t) to estimate l/Ti. The param- 
eters are set as follows: the onsite interaction Vu — —2t, 
the chemical potential 11 = —t, and the system size 
Lj, X Lj, = 64 X 64. The delta function 6{x) is approx- 
imated by 6(x) — {1 / T:)rj / {x^ -f rf ) with the smearing 
factor 7y = O.Oli. We adopt the shifted BiCG method as 
a iterative linear solver. As shown in Fig.[TTl the nuclear 
magnetic relaxation rate can be successfully reproduced 
by the SS method. We note that the discrete energy 
levels due to the finite size system cause the relatively 
smaller Hebel-Slichter peak below Tc- 



F. Computational costs in the eigenvalue problem 
with the SS-method in a vortex lattice system 

Now, let us evaluate the computational costs of the 
SS method. We measure the elapsed time from read- 
ing the Hamiltonian matrix constructed by the poly- 
nomial expansion scheme to finishing the SS-method 
in a Lj; X Ly square lattice s-wave superconductor at 
T = 0.04t {Lx = Ly). We use the contour F with 7 = 
and p = O.li, and the physical parameters are the same 
as in Fig. [S] For the measurement, we use a desktop 
computer with only one CPU core (Intel Xeon X5550 
2.66GHz). As shown in Fig. [T^l the elapsed time of the 
SS method grows in an 0{N) manner with increasing 
the system size N = 2{Lx x Ly). The computational 
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System size N 



FIG. 12. (Color online) System-size dependence of the 
elapsed time with one CPU core for obtaining the eigenval- 
ues {—O.lt < < O.lt) with the SS-method, for an s-wave 
superconductor at T = 0.04t, on a x Ly square lattice 
(Lx = Ly). The system size is A'^ = 2Lx x Ly. 



costs are roughly estimated by 0{msN) In all the calcula- 
tion of this subsection, the energy domain is restricted to 
—O.lt < ei < Q.lt. As a result, the number of the eigen- 
values iris is independent of the spatial size N . Indeed, 
we find that only the bound states in vortices are rele- 
vant to this narrow energy window. If an energy domain 
is wide (large p), rus increases. In this case, the computa- 
tional cost predominantly depends on the Gram- Schmidt 
orthonormalization procedure to construct Q. Then, we 
find that the cost is estimated by 0{Nm'^). When one 
tries to obtain all the eigenvalues (i.e., to^ = N) with a 
wide energy domain, the cost of the SS method becomes 
0{N^). This result is equivalent to the one in the full 
diagonalization method. However, the wide energy do- 
main is easily divided into small energy domains. For 
example, with using TVs domains with Ng parallel com- 
putation with Ns CPU cores, the elapsed time reduces to 
1/A^s- This is a big advantage of the SS method. 



VI. CONCLUSION 

We proposed the fast efficient method on the basis 
of the SS-method and the polynomial expansion scheme 
to calculate the eigen-pairs and the dynamical correla- 
tion functions in the BdG scheme of superconductivity. 
The polynomial expansion scheme enables us to solve 
the gap-equations self consistently, with large scale par- 
allel computations. With the use of the SS method, 
one can solve the problem of finding several eigenvalues 
in a given domain and their corresponding eigenvectors. 
We applied the present approach to the calculations of 
various physical quantities including the eigenvalues dis- 



tribution of the BdG Hamiltonian in a vortex lattice, 
magnetic-field-dependence of the thermal conductivity, 
and temperature-dependence of the nuclear magnetic re- 
laxation rate. We stress that most of the calculations 
were performed, changing the system size. This is quite 
important for developing a theoretical tool to predict 
physical behaviors in nano-scale superconductors from a 
microscopic theory. 
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Appendix A: Iterative refinement 

In Sec. IIV CI we remark that larger subspace makes 
more higher accuracy but needs more heavier computa- 
tional costs. A problem may occur when we choose a 
small value of L to avoid the use of a large subspace. 
If some residuals of the obtained approximate eigenval- 
ues and eigenvectors are not small enough for a given 
tolerance, we can brush up the resulting approximate 
eigenvalues and eigenvectors. We propose two ways to 
brush up the eigen pairs with the choice of the appro- 
priate source matrix V dominantly constructed by the 
information in a given domain F. 

First, one can use the source matrix V as 

V - 5o, (Al) 
Then, we implement r iterations of Pt{A) on V ^ 

Si"-^ = PriA)st'\ si^^^V. (A2) 

Using 5*0 , we construct a refined matrix including 

a higher moment vector, s'^'' —A'^Pr{A)SQ^ ^\ In a 
simulation, numerical quadrature is used for calculat- 
ing PriA) and A'=Pr(^), as seen in Sec. HVB] Thus, 
performing the singular-value decomposition of S^^^ — 

^ (r) " (r) 

{Sq , . . . , SIj[_^}, we evaluate a refined effective rank mg. 

Second, we can brush up the resulting approximate 
eigenvalues and eigenvectors by setting the source matrix 
as 

V = {xi,--- ,x,nJC, (A3) 

where C £ C™"^^ whose elements are random numbers 
in (—1,1), and Xi, - ■ ■ , are the selected eigenvectors 
that are regarded as the approximate eigenvectors with 
respect to the eigenvalues inside F. Using this V, we 
reevaluate S and mg in Sec. IIVBI 
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Appendix B: Estimation of the trace 

We show that the trace of an n x n matrix A can be 
estimated by 



(Bl) 



fc=i 



with random vectors with entries ±1. If the vectors Vk 
have entries ±1, the right-hand side in the above equation 
is expressed as 



S ^ 71 S 



(B2) 



On the average, the coefficient of Aij in the above expan- 
sion will converge to zero provided that the components 
of the vectors v'' have balanced ± signs. 



Appendix C: Thermal conductivity 

We derive the expression of the thermal conductivity 
in terms of the solution of the BdG equation. Let us con- 
sider the case when a temperature gradient exists along 
X-axis on a two-dimensional x Ly lattice with lattice 
constant a. The electronic thermal conductivity per vol- 
ume of a superconductor associated with heat flux along 
X-axis is given in Eq. (|27|) . The current-current correla- 
tion function with respect to the energy flux is 

p,,(r,o) = (r,[Ju,x(r)Ju,.(o)]) ^^Y. e-'''"'^PU^n„,) 

(CI) 

The Heisenberg operator of energy flux— along x-axis is 
'dc]. 



n -r — D* r 



t dCj^a 



(C2) 



with Dx,ij = Si+i^j{—iat)e^ ^^ and Ix = {1,0) = a/a. 
The link variable e represents the contribution of the 
magnetic field, 6ij = {■7T/(j)o)a ■ A[{ri -\- rj)/2], with the 
flux quantum 0o and the vector potential A. The matrix 
Dx{= (Dx^ij)) is related to the momentum operator on a 
square lattice—. 

This current-current correlation function may be 
rewritten as the form of a two-particle Green's func- 
tion. 



28 



X {Tr [^1 (rO^fc (ri (^^ (rs)] ) , (C3) 

with t{ — >• Ti -I- and — > T2 -I- 0. We neglected 
some terms associated with the action of the 
imaginary-time derivative on the time ordering oper- 
ator, according to the discussion by Ambegaokar and 
Teword^- . Here, tp and -0^ are defined as, respectively, 
= (ct-,C4_)'^ and V'^ = (c^, cj), with c^r = (ci^a, ■ ■ 

and Ca = (c| ^, .. .)"'". This convention corresponds to 
the Nambu representation. The indices a, b, c, and d run 
from 1 to 2LxLy. The 2LxLy x 2LxLy matrix J^^^ is 
defined as 











(C4) 



with the LxLy x LxLy matrix 

J^''^^r',^r) = -i{dr'Dx " drOl). Thc formula- 
tion in Secini which is developed in a more general 
representation for a fermion mean-field theory, is 
straightforwardly rewritten in terms of the Nambu 
representation. 

Now. let us derive the expression of the thermal 
conductivity. First, we evaluate the imaginary part of 
Pxx{i^m ^ + iO). Next, we expand the resultant for- 
mula up to fi, to take the limit 57 — 0. We note that 
lm[Pxx{^ = 0)] = 0. Then, we obtain Eq. (gT]). The 
2LxLy X 2LxLy matrix V'-^^ in Eq. (l27l) is defined as 



{Dx 







(C5) 
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